# -------------------------------------#
#         FOX NEWS & COVID-19          #
#         SURVEY ANALYSIS              #
# -------------------------------------#

## Clear environment and return memory to the operating system
rm(list = ls())
gc()

## Installing and loading the needed packages
packages_required <- c('stargazer', 'foreign', 'lmtest', 'haven', 'expss', 'dplyr',
                       'jtools', 'ggstance', 'coefplot',  'ggplot2', 'texreg', 'broom', 'patchwork',
                       'sandwich', 'cowplot', 'tidyverse', 'Ecdat', 'estimatr', 'usethis', 'devtools')
for (package in packages_required) {
  if (!(package %in% installed.packages())) {
    install.packages(package)
  }
}
Packages <- c('stargazer', 'foreign', 'lmtest', 'haven', 'expss', 'dplyr', 'patchwork',
              'jtools', 'ggstance', 'coefplot',  'ggplot2', 'texreg', 'broom',
              'sandwich', 'cowplot', 'tidyverse', 'Ecdat', 'estimatr', 'usethis', 'devtools')
lapply(Packages, library, character.only = TRUE)  

# Dataset
dataR3 <- read_dta("")  # Read in the dta file generated with the rowdata dofile.
datat<-as.tibble(dataR3) #a tibble version
glimpse(datat)

#A subset of the dataframe - only respondents who pass the screener
datat <- subset(dataR3, pass_screeners==1)
glimpse(datat)

#-----------------------------------------------------------------------------------------------#
#      TABLE 1 --  The effect of watching Fox News, CNN & MSNBC on three dependent variables    #
#                  including state fixed effects and clustering by state - subset of cabletv    # 
#-----------------------------------------------------------------------------------------------#

#subset of the dataframe - only responsents who watch cable tv chanels on a daily or weekly basis
glimpse(datat)
table(datat$cabletvdum)
levels(as_factor(datat$cabletv))
cabletv_data <- subset(datat, cabletvdum==1)
glimpse(cabletv_data)

#outcome 1: changing behavior - whether a respondent changes his behavior in the last week due to the coronavirus
matrixbhv <- lm_robust(matrixbhv_mean~foxnews_watch+cnn_watch+msnbc_watch+
                         age+income+female+education+factor(race)+factor(PID3),
                       data=datat,clusters=State,
                       fixed_effects=~State)
summary(matrixbhv)

#without PID
matrixbhv_whithoutpid <- lm_robust(matrixbhv_mean~foxnews_watch+cnn_watch+msnbc_watch+
                                     age+income+female+education+factor(race),
                                   data=datat,clusters=State,
                                   fixed_effects=~State)
summary(matrixbhv_whithoutpid)

#outcome 2: economic-health tradeoff - whether a respondent prefers investing in public health over the economy in managing the Corona outbreak
health_tradeoff <- lm_robust(ecohealth_tradeoff2~foxnews_watch+cnn_watch+msnbc_watch+
                               age+income+female+education+factor(race)+factor(PID3),
                             data=datat,clusters=State,
                             fixed_effects=~State)
summary(health_tradeoff)
#without PID
health_tradeoff_whithoutpid <- lm_robust(ecohealth_tradeoff2~foxnews_watch+cnn_watch+msnbc_watch+
                                           age+income+female+education+factor(race),
                                         data=datat,clusters=State,
                                         fixed_effects=~State)
summary(health_tradeoff_whithoutpid)

#outcome 3: Hydroxychloroquine - whether a respondent believes that the drug HC is an effective treatment for the coronavirus 
Hydroxychloroquine <- lm_robust(Hydroxychloroquine~foxnews_watch+cnn_watch+msnbc_watch+
                                  age+income+female+education+factor(race)+factor(PID3),
                                data=datat,clusters=State,
                                fixed_effects=~State)
summary(Hydroxychloroquine)
#without PID
Hydroxychloroquine_whithoutpid  <- lm_robust(Hydroxychloroquine~foxnews_watch+cnn_watch+msnbc_watch+
                                               age+income+female+education+factor(race),
                                             data=datat,clusters=State,
                                             fixed_effects=~State)
summary(Hydroxychloroquine_whithoutpid)

## define the 3 regression models as a dataframe
Model1 <-  tidy(health_tradeoff, conf.int = TRUE)
Model2 <-  tidy(Hydroxychloroquine, conf.int = TRUE)
Model3 <-  tidy(matrixbhv, conf.int = TRUE)
allmodels <- bind_rows(Model1, Model2, Model3, .id="Outcome")
glimpse(allmodels)

allmodels$Outcome <- factor(allmodels$Outcome, labels = c("Health-Economic Tradeoff",
                                                          "Hydroxychloroquine is Effective",
                                                          "Changing Behavior"))
allmodels <- subset(allmodels, term=="foxnews_watch" | 
                      term=="cnn_watch" | 
                      term=="msnbc_watch")

allmodels$Newschannel <- factor(allmodels$term, labels = c("CNN",
                                                           "Fox News",
                                                           "MSNBC"))

allmodels$Newschannel2 <- factor(allmodels$Newschannel, levels = c("Fox News", 
                                                                   "CNN", 
                                                                   "MSNBC"))

#------------------------------------------------------------#
#                        COEFPLOT 1                          # 
#------------------------------------------------------------#
fullsample<-
  ggplot(allmodels, 
         aes(factor(Outcome), estimate, 
             color=factor(Newschannel2), 
             group=factor(Newschannel2), 
             shape=factor(Newschannel2)))+ 
  scale_colour_manual(values = c("#175fad", '#bf1d00', "#000000"))+
  scale_y_continuous(breaks = c(-0.2, -0.1, 0.0, 0.1, 0.2), limits = c(-0.3, 0.3))+
  geom_hline(yintercept=0, lty=2, lwd=0.7, size=0.05, colour = "grey") +
  geom_errorbar(stat = "identity", alpha = 0.5, 
                position = position_dodge(width = 0.15),
                aes(ymin=estimate - 1.96*std.error, ymax=estimate + 1.96*std.error), 
                lwd=1, width=0) +
  geom_errorbar(stat = "identity", alpha = 0.8, 
                position = position_dodge(width = 0.15),
                aes(ymin=estimate - 1.64*std.error, ymax=estimate + 1.64*std.error), 
                lwd=1, width=0, size=3) +
  geom_point(stat = "identity", 
             alpha = 1, 
             position = position_dodge(width = 0.15), 
             size = 2) +
  ggtitle("") +
  coord_flip() +
  labs(x = "", y = "", title = "", subtitle = "", 
       color="News Channel", face = "bold", 
       shape="News Channel", face = "bold") +
  theme_minimal() +
  theme(legend.position = "right")+
  theme(plot.title = element_text(size = 10, face = "bold"),
        axis.text.y = element_text(face = "bold"),
        axis.text.x = element_text(face = "bold"))
fullsample

#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#
#      TABLE 2 --  Fox News Programs & Timing of COVID-19 Related Behaviors                                                                 # 
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#

#controls -- complete sample
glimpse(datat) #only pass screeners 
model_fullsample <- lm_robust(matrixbhv_mean~Hannity_watch+FoxFriends_watch+Tucker_watch+
                                age+income+female+education+factor(race),
                              data=datat,
                              clusters=State,
                              fixed_effects=~State)
summary(model_fullsample)

#controls -- Republicans only
data_Rep <- subset(datat, PID3==1)
model_Rep <- lm_robust(matrixbhv_mean~Hannity_watch+FoxFriends_watch+Tucker_watch+
                         age+income+female+education+factor(race),
                       data=data_Rep,
                       clusters=State,
                       fixed_effects=~State)
summary(model_Rep)

#controls-  Republicans and foxnews whatcher
data_Rep_foxnews <- subset(data_Rep, foxnews_watch==1)
model_Rep_foxnews <- lm_robust(matrixbhv_mean~Hannity_watch+FoxFriends_watch+Tucker_watch+
                                 age+income+female+education+factor(race),
                               data=data_Rep_foxnews,
                               fixed_effects=~State)
summary(model_Rep_foxnews)
#controls-  Republicans and conservative
data_Repcons <- subset(data_Rep, libcons3==3)
model_Repcons <- lm_robust(matrixbhv_mean~Hannity_watch+FoxFriends_watch+Tucker_watch+
                             age+income+female+education+factor(race),
                           data=data_Repcons,
                           fixed_effects=~State)
summary(model_Repcons)


#----------------------------------------#
#                 PLOT 2                 #
#----------------------------------------#

model_fullsample <-  tidy(model_fullsample, conf.int = TRUE)
model_Rep <-  tidy(model_Rep, conf.int = TRUE)
model_Rep_foxnews <-  tidy(model_Rep_foxnews, conf.int = TRUE)
model_Repcons <-  tidy(model_Repcons, conf.int = TRUE)

library(dplyr)
allModelFrame2 <- bind_rows(model_fullsample, model_Rep, model_Rep_foxnews, model_Repcons, .id="model")
glimpse(allModelFrame2)
allModel2 <- subset(allModelFrame2, term=="FoxFriends_watch"|
                      term=="Hannity_watch"|
                      term=="Tucker_watch")

glimpse(allModel2)
allModel2$model <- factor(allModel2$model,
                          levels = c(1,2,3,4),
                          labels = c("Full Sample",
                                     "Republicans Only",
                                     "Republicans who are Fox News Viewers",
                                     "Republicans who are Conservatives"))
allModel2$term2 <- factor(allModel2$term,
                          labels = c("Fox & Friends",
                                     "Hannity",
                                     "Tucker Carlson Tonight"))
allModel2$term3[allModel2$term2 == "Hannity"] = 3
allModel2$term3[allModel2$term2 == "Fox & Friends"] = 1
allModel2$term3[allModel2$term2 == "Tucker Carlson Tonight"] = 2
allModel2$term3 <- factor(allModel2$term3,
                          levels=c(1,2,3),
                          labels = c("Fox & Friends",
                                     "Tucker Carlson Tonight",
                                     "Hannity"))
glimpse(allModel2)

# Specify the width of your confidence intervals
interval1 <- -qnorm((1-0.9)/2)  # 90% multiplier
interval2 <- -qnorm((1-0.95)/2)  # 95% multiplier


program_plot1<-
  ggplot(allModel2, 
         aes(factor(term3), estimate, 
             color=factor(term3), 
             group=factor(term3), 
             shape=factor(term3)))+ 
  scale_colour_manual(values = c("#175fad", '#bf1d00', "#000000"))+
  scale_y_continuous(breaks = c(-0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8))+
  geom_hline(yintercept=0, lty=2, lwd=0.7, size=0.05, colour = "grey") +
  geom_errorbar(aes(ymin=estimate - 1.96*std.error, ymax=estimate + 1.96*std.error), 
                lwd=0.9, width=.0, size=0.35, alpha = 0.8) +
  geom_errorbar(aes(ymin=estimate - 1.645*std.error, ymax=estimate + 1.645*std.error), 
                lwd=1.5, width=.0, size=0.35, alpha = 0.5) +
  geom_point(size = 3.3) +
  ggtitle("") +
  coord_flip() +
  labs(x = "", y = "", title = "", subtitle = "") +
  facet_wrap(~factor(model))+
  theme_minimal() +
  theme(legend.position = "none",
        panel.spacing = unit(2, "lines")) +
  theme(plot.title = element_text(size = 10, face = "bold")) +
  theme(strip.text = element_text(face = "bold"))
program_plot1


#----------------------------------------------------------------------#
#      TABLE 2 --  Fox News Programs & health_tradeoff                 # 
#----------------------------------------------------------------------#

#controls -- complete sample
model_fullsample2 <- lm_robust(ecohealth_tradeoff2~Hannity_watch+FoxFriends_watch+Tucker_watch+
                                 age+income+female+education+factor(race),
                               data=datat,
                               clusters=State,
                               fixed_effects=~State)
summary(model_fullsample2)

#controls -- only Republicans
data_Rep <- subset(datat, PID3==1)
model_Rep2 <- lm_robust(ecohealth_tradeoff2~Hannity_watch+FoxFriends_watch+Tucker_watch+
                          age+income+female+education+factor(race),
                        data=data_Rep,
                        clusters=State,
                        fixed_effects=~State)
summary(model_Rep2)

#controls-  Republicans and foxnews whatcher
data_Rep_foxnews <- subset(data_Rep, foxnews_watch==1)
model_Rep_foxnews2 <- lm_robust(ecohealth_tradeoff2~Hannity_watch+FoxFriends_watch+Tucker_watch+
                                  age+income+female+education+factor(race),
                                data=data_Rep_foxnews,
                                fixed_effects=~State)
summary(model_Rep_foxnews2)

#controls-  Republicans and conservative
data_Repcons <- subset(data_Rep, libcons3==3)
model_Repcons2 <- lm_robust(ecohealth_tradeoff2~Hannity_watch+FoxFriends_watch+Tucker_watch+
                              age+income+female+education+factor(race),
                            data=data_Repcons,
                            fixed_effects=~State)
summary(model_Repcons2)


#----------------------------------------------------------------------#
#      TABLE 2 --  Fox News Programs & Hydroxychloroquine              # 
#----------------------------------------------------------------------#

#controls -- complete sample
model_fullsample3 <- lm_robust(Hydroxychloroquine~Hannity_watch+FoxFriends_watch+Tucker_watch+
                                 age+income+female+education+factor(race),
                               data=datat,
                               clusters=State,
                               fixed_effects=~State)
summary(model_fullsample3)

#controls -- only Republicans
data_Rep <- subset(datat, PID3==1)
model_Rep3 <- lm_robust(Hydroxychloroquine~Hannity_watch+FoxFriends_watch+Tucker_watch+
                          age+income+female+education+factor(race),
                        data=data_Rep,
                        clusters=State,
                        fixed_effects=~State)
summary(model_Rep3)

#controls-  Republicans and foxnews whatcher
data_Rep_foxnews <- subset(data_Rep, foxnews_watch==1)
model_Rep_foxnews3 <- lm_robust(Hydroxychloroquine~Hannity_watch+FoxFriends_watch+Tucker_watch+
                                  age+income+female+education+factor(race),
                                data=data_Rep_foxnews,
                                fixed_effects=~State)
summary(model_Rep_foxnews3)

#controls-  Republicans and conservative
data_Repcons <- subset(data_Rep, libcons3==3)
model_Repcons3 <- lm_robust(Hydroxychloroquine~Hannity_watch+FoxFriends_watch+Tucker_watch+
                              age+income+female+education+factor(race),
                            data=data_Repcons,
                            fixed_effects=~State)
summary(model_Repcons3)

